﻿/******************************************************************************
 * 
 * Announce: CSharpKit, Basic algorithms, components and definitions.
 *           Copyright (C) ShenYongchen.
 *           All rights reserved.
 *   Author: 申永辰.郑州 (shenyczz@163.com)
 *  WebSite: http://github.com/shenyczz/CSharpKit
 *
 * THIS CODE IS LICENSED UNDER THE MIT LICENSE (MIT).
 * THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF 
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
 * IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
 * PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
 * 
******************************************************************************/

using System;
using System.Numerics;

namespace CSharpKit.Numerics.LinearAlgebra.Factorization
{
    /// <summary>
    /// SVD - 奇异值分解（singular value decomposition）
    /// <a href="http://en.wikipedia.org/wiki/Binomial_distribution">Wikipedia - Singular value decomposition</a>.
    /// </summary>
    /// <remarks>
    /// 理论描述：<para/>
    /// 
    /// Suppose M is an (m x n) matrix whose entries are real numbers.
    /// Then there exists a factorization of the form M = UΣVT where: <br/>
    /// - U  is m-by-m unitary matrix（酉阵） <br/>
    /// - Σ is m-by-n diagonal matrix with nonnegative real numbers on the diagonal; <br/>
    /// - VT denotes transpose of V, an n-by-n unitary matrix; <para/>
    /// 
    /// The Σ(i,i) on the diagonal is the SVD of matrix M.<br/>
    /// A common convention(惯例) is to order the diagonal entries Σ(i,i) in descending order.
    /// In this case, the diagonal matrix Σ is uniquely determined
    /// by M (though the matrices U and V are not). The diagonal entries of Σ are known as the singular values of M.
    /// </remarks>
    /// <typeparam name="T">Supported data types are double, single, <see cref="Complex"/>, and <see cref="Complex32"/>.</typeparam>
    public abstract class Svd<T> : ISolver<T>
        where T : struct, IFormattable, IEquatable<T>
    {
        protected Svd(Vector<T> s, Matrix<T> u, Matrix<T> vt, bool vectorsComputed)
        {
            S = s;
            U = u;
            VT = vt;

            VectorsComputed = vectorsComputed;

            _lazyW = new Lazy<Matrix<T>>(ComputeW);
        }





        readonly Lazy<Matrix<T>> _lazyW;

        /// <summary>
        /// 指示在SVD分解过程中是否计算了U和VT矩阵。
        /// Indicating whether U and VT matrices have been computed during SVD factorization.
        /// </summary>
        protected readonly bool VectorsComputed;


        /// <summary>
        /// Gets the singular values (Σ) of matrix in ascending value.
        /// </summary>
        public Vector<T> S { get; private set; }

        /// <summary>
        /// Gets the left singular vectors (U - m-by-m unitary matrix)
        /// </summary>
        public Matrix<T> U { get; private set; }

        /// <summary>
        /// Gets the transpose right singular vectors (transpose of V, an n-by-n unitary matrix)
        /// </summary>
        public Matrix<T> VT { get; private set; }


        /// <summary>
        /// the singular values as a diagonal <see cref="Matrix{T}"/>.
        /// </summary>
        public Matrix<T> W => _lazyW.Value;


        public abstract T Determinant { get; }
        public abstract int Rank { get; }
        /// <summary>
        /// Gets the two norm of the <see cref="Matrix{T}"/>.
        /// </summary>
        public abstract double L2Norm { get; }
        /// <summary>
        /// Gets the condition number <b>max(S) / min(S)</b>
        /// </summary>
        /// <returns>The condition number.</returns>
        public abstract T ConditionNumber { get; }

        private Matrix<T> ComputeW()
        {
            var rows = U.RowCount;
            var columns = VT.ColumnCount;
            var result = Matrix<T>.Builder.SameAs(U, rows, columns);
            for (var i = 0; i < rows; i++)
            {
                for (var j = 0; j < columns; j++)
                {
                    if (i == j)
                    {
                        result.At(i, i, S[i]);
                    }
                }
            }

            return result;
        }



        #region ISolver<T>

        public Matrix<T> Solve(Matrix<T> input)
        {
            if (!VectorsComputed)
            {
                throw new InvalidOperationException(Resources.Numerics.SingularVectorsNotComputed);
            }

            var x = Matrix<T>.Builder.SameAs(U, VT.ColumnCount, input.ColumnCount, fullyMutable: true);
            Solve(input, x);
            return x;
        }
        public abstract void Solve(Matrix<T> input, Matrix<T> result);

        public Vector<T> Solve(Vector<T> input)
        {
            if (!VectorsComputed)
            {
                throw new InvalidOperationException(Resources.Numerics.SingularVectorsNotComputed);
            }

            var x = Vector<T>.Builder.SameAs(U, VT.ColumnCount);
            Solve(input, x);
            return x;
        }
        public abstract void Solve(Vector<T> input, Vector<T> result);

        #endregion


        //}}@@@
    }


}

